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We study the effects of anisotropic hopping amplitudes on quantum phases of ultracold fermions 
in optical lattices described by the repulsive Fermi-Hubbard model. In particular, using dynamical 
mean-field theory (DMFT) we investigate the dimensional crossover between the isotropic square 
and the isotropic cubic lattice. We analyze the phase transition from the antiferromagnetic to the 
paramagnetic state and observe a significant change in the critical temperature: Depending on the 
interaction strength, the anisotropy can lead to both a suppression or increase. We also investigate 
the localization properties of the system, such as the compressibility and double occupancy. Using 
the local density approximation in combination with DMFT we conclude that density profiles can 


be used to detect the mentioned anisotropy-driven 
PACS numbers: 71.10.Fd, 75.50.Ee, 67.85.-d, 37.10.Jk 

I. INTRODUCTION 

Ultracold atoms in optical lattices offer not only a clean 
and controlled quantum simulator for electronic solid- 
state materials, but also constitute interesting many- 
body systems in their own right, exhibiting a rich variety 
of physical phenomena. The high degree of tunability 
of the parameters allows to study the system in various 
regimes and thus to make direct quantitative compar¬ 
isons with the predictions of theoretical models describ¬ 
ing quantum many-body phenomena, such as the promi¬ 
nent Hubbard model pQ . A periodic potential is an essen¬ 
tial ingredient of these models because it largely deter¬ 
mines the electronic behavior in solid materials. Thus, 
the implementation of optical lattices in cold atom se¬ 
tups Eiia set a milestone in this experimental field. Since 
then, a substantial progress was made towards realizing 
the idea of a universal quantum simulator by preparing 
two-component mixtures of ultracold fermionic atoms in 
optical lattices. In particular, a fermionic Mott insu¬ 
lator was realized and detected BUS, short-range anti¬ 
ferromagnetic correlations were observed and effectively 
measured 00, and recently, the single-site resolution of 
fermions in optical lattices was experimentally attained 

mm- 

The lattice geometry is a crucial characteristic of these 
systems, since it essentially affects all the other physical 
properties. Experimentally, the geometry of an optical 
lattice is determined, among other things, by the spa¬ 
tial arrangement of the laser beams, and the possibil¬ 
ity for the trapped atoms to tunnel within the lattice is 
then tuned via the laser intensity mm- In particular, 
starting from an isotropic cubic lattice and gradually in¬ 
creasing the laser intensity along the ^-direction results 
in a system consisting of separated layers of square lat¬ 
tices parallel to the xy-plane and thus effectively corre¬ 
sponds to a dimensional crossover from 3d to 2d (another 
type, from 3d to Id, was studied in the related context in 


transitions. 


Refs. 0 |T3] ) • We analyze the effects of this dimensional 
crossover for a fermionic system with repulsive interac¬ 
tions. In particular, we study the magnetically ordered 
phases of the system and its thermodynamic properties. 
Employing a combination of DMFT and local density ap¬ 
proximation (LDA+DMFT), we perform calculations for 
the case of an additional external potential (harmonic 
trap) and analyze the real-space particle distribution as 
well as the entropy of the system in different parameter 
regimes. 


II. SYSTEM AND MODEL 

The theoretical description is given by the Fermi- 
Hubbard model and the anisotropy is constituted in the 
ratio of the hopping parameters t z /t, where the hopping 
parameters in the xy-plane are set equal (t x = t y = t) 
and t z varies between zero and t (0 < t z /t < 1). We 
consider a Fermi-Hubbard Hamiltonian of the following 
type: 

H =-'^2'^2Uj(c\ a c jcr +h.c.) + U Y1 

(ij) a i 

i a 

where the notation (ij) indicates a summation over 
nearest-neighbor sites, tij is the hopping amplitude of 
the fermions with values for each spatial direction spec¬ 
ified above and U is the magnitude of the on-site re¬ 
pulsive (U > 0) interaction of the two different species 
(or hyperfine states) a = The translational 

eigenstates in the periodic lattice potential Vi at ( r) = 
Vq°^ cos 2 (ha) with lattice depth Vq“^ can be rep- 

cx.=x,y,z 

resented in terms of Wannier orbitals. With i/j (r — r/ 
being a single Wannier function localized at site i the 
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parameters for hopping and onsite interaction can be ex¬ 
pressed in terms of overlap integrals: 


Uj = - J d 3 r ip*(r - I-*) V>(r - r,), 


U = 


Att h 2 a s 
m 


1 


d 3 r |^(r)| 4 . 


The Weiss function can be expressed as Go^{ioJ n ) = 
iu> n + fi — A a (ico n ), with the fermionic Matsubara fre¬ 
quency u) n = (2 n +1 )it//3, where n is an integer number 
and /3 is the inverse temperature, and the hybridization 
function Awhich is determined by the Anderson 
parameters, 


A a (i(jJ n ') 
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In expression (JT|) , cj CT (ci CT ) is the corresponding creation 
(annihilation) operator of species with (pseudo)spin er, 
with corresponding densities and h q. (= cj^c^). 

Vi is the external (e.g., harmonic) potential at lattice site 
i, and p is the chemical potential of the atoms. Note that 
we set the harmonic potential and the optical lattice to be 
independent of the atomic species. The Hamiltonian 0 
implies a single-band approximation; in other words, we 
consider the case of a sufficiently strong lattice potential, 
Vq > 5 E r , with E r = H 2 k 2 /(2m) being the recoil energy 
of the fermions. 


III. METHOD 

Our analysis is based on the dynamical mean-field the¬ 
ory (DMFT). In this approach the system is reduced to 
a local many-body problem described as a single lattice 
site (“impurity”) coupled to an “external bath”, i.e. the 
spatial fluctuations are frozen, but the local dynamics of 
the system is fully preserved. The self-consistency con¬ 
ditions of DMFT can be derived in various ways. The 
approach introduced by A. Georges and G. Kotliar p~4) 
uses the Anderson impurity model (15], considering the 
single site as an “impurity orbital” and the bath as a 
“conduction band”: 

Id am ^ ^ T ^ ) V[fj (d^C/j A c^.h/ CT ) 

l,a l,a 

+C/nfnf-/r(nf+ fif) (2) 

is the corresponding Hamiltonian with the fermionic cre¬ 
ation and annihilation operators cj., c a describing the lo¬ 
cal degrees of freedom and a/ lcr , di a describing a set of non¬ 
interacting fermions representing the degrees of freedom 
of the effective bath acting on the site i. Hence, the first 
term in the Hamiltonian (|2j) describes the effective bath, 
the last two terms describe the impurity site, and the sec¬ 
ond term constitutes the coupling. Bath orbitals are de¬ 
noted by the index l, ei„ and Vi a are the so-called Ander¬ 
son parameters. Conveniently, the DMFT equations are 
expressed in the functional integral formalism. The key 
quantities in this notation are the so-called Weiss func¬ 
tion Go,a(iu n ) and the local propagator G a {iui n ), which 
are related to the self-energy Yi a (iuj n ) by the Dyson equa¬ 
tion: 


(3) 


The self-consistency condition relating the dynamical 
mean-field A a (iuj n ) to the local Green’s function G a (ioj n ) 
is given by 


G tjiiuin) 
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Impurity solver. The Anderson impurity model can 
be solved by various techniques. We choose the exact 
diagonalization (ED) impurity solver [TBj. Within this 
solver the number of orbitals in the effective Anderson 
model is considered to be finite. The algorithm essen¬ 
tially consists of three steps: (i) The Weiss function is 
approximated by a discretized version. This replacement 
can be considered as a projection on a restricted func¬ 
tional subspace, which is the most involved part of the 
diagonalization procedure. The existing methods for the 
implementation of this projection are described in refer¬ 
ence El- ( ii) The Anderson Hamiltonian ([2| restricted 
to a finite number ( n s ) of orbitals is diagonalized ex¬ 
actly and the local propagator G a (iuj n ) is computed, (iii) 
The self-consistency condition yields a new Weiss func¬ 
tion Go,a, which in turn is approximated by a function 
Gq 1 ^ with a new set of Anderson parameters ii a and Vj CT . 
The process is iterated until a converged set of param¬ 
eters is reached. The full numerical diagonalization is 
currently feasible up to n s = 7. Our results are based on 
calculations with n s = 5. 

Two-sublattice DMFT. Bipartite structures, such as 
states with antiferromagnetic order, can be described 
in an appropriate way by introducing two sublattices. 
Within the DMFT approach one then needs to solve the 
impurity problem twice on two adjacent sites of the orig¬ 
inal lattice. The self-consistency condition for this case 
is given by 


G;(^) = c/±gAW (?) 

with = iui n A /i — , and the sublattice indices a = 

A,B and its opposite a = B,A. In the balanced case 
(i.e. no imbalance in the hopping amplitudes or chemical 
potentials between two spin components El), there is 
an additional symmetry present, () 4 = , which allows 

to reduce the computations related to the solution of the 
impurity problem to a single site. 


Yti a (iiVn) — @0,(r(^ n ') dd u (ioJ n ). 
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FIG. 1. (Color online) Non-interacting density of states pro¬ 
files for t x /t = ty/t = 1 and different G-values. 


Note that (apart from the bipartite structure of the 
self-consistency condition) the geometry of the lattice 
enters our numerical studies only through the non¬ 
interacting density of states D(e). As can be seen in 
Fig .0 this quantity depends strongly on the possible 
anisotropies in the system. In Appendix [A] we derive an¬ 
alytical expressions for D(e). The evaluation of elliptic 
integrals involved in these expressions was done by the 
AGM method and is presented in Appendix |ES[ 


(a) 
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FIG. 2. (Color online) (a) Phase diagrams for the Neel- 
ordered phases in lattices with the hopping parameter t z /t = 
0, 0.5, 1 obtained by DMFT at half filling (fj = U/ 2). The 
antiferromagnetic order parameter, staggered magnetization 
m = |n;-|- — is color-coded in the T/U- plane, (b) Depen¬ 
dence of the critical temperature T c on the hopping parameter 
t z for different interaction strengths U. The critical tempera¬ 
ture is nearly constant for U = 6t and it increases (decreases) 
for larger (smaller) values of U. 


IV. RESULTS 

In this section we discuss how the introduced spatial 
anisotropy in hopping amplitudes influences the main 
many-body properties of the system. 

Preceding the following discussion it should be noted 
that DMFT calculations for a 2d system predict anti¬ 
ferromagnetic long-range order (LRO), whereas it is well 
known that according to the Mermin-Wagner-Hohenberg 
theorem LRO can arise in the isotropic (SU(2) sym¬ 
metric) Hubbard model at finite temperatures only for 
d > 2 [18]. Thus, in low-dimensional systems (d < 2) 
LRO for this model is only possible at T = 0. However, 
the DMFT results for these systems do have a physical 
relevance - the calculated Neel temperature gives a good 
quantitative estimate for the boundary at which (short- 
ranged) antiferromagnetic correlations arise in the sys¬ 
tem TT)|. 

A. Magnetically ordered phases 

As expected on the basis of known limiting cases 
for isotropic hopping in square (t z /t = 0) and cubic 


(t z /t = 1) lattice geometries, we observe a continuous 
change of the ordered phase region with an increase of t z 
as shown in Fig. [2ja) . In particular, at strong coupling 
(U > 5 1) the critical temperature increases significantly, 
in the region of the intermediate interaction strength U/t 
the phase transition line stays approximately unchanged, 
but at low U/t the antiferromagnetic phase is more sup¬ 
pressed. To demonstrate these effects we plot the criti¬ 
cal temperature as a function of t z for different interac¬ 
tion strengths (Fig. [2jb)) . With increasing values of t z , 
the critical temperature decreases for weak interactions, 
rises for strong interactions and stays approximately un¬ 
changed for U/t = 6. This behavior becomes plausible 
when we consider the known analytical results for the 
limiting cases. 

In the case U <C Zt the critical temperature is given 
by 

T c oc VZte~ cV v 

with the coordination number Z and a positive constant 
c [20] . In this regime, the exponential suppression dom¬ 
inates over the linear growth with increase of the coor¬ 
dination number Z which corresponds to increase of t z 
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Interaction U/t 

FIG. 3. (Color online) Contour lines (/t e = 0.01 /t) (solid) 
and magnetic transition lines (dashed) for t z /t = 0,0.5,1. 
For larger t z the crossover to the insulating regime is shifted 
towards larger interaction strength values. 


in our case. For the opposite case of strong coupling 
U Zt the Heisenberg spin model is applicable. Within 
the mean-field description T c is proportional to the mag¬ 
netic coupling constant J, 



As one can see from Fig. [2j)b) , the mean-field relation 
between the critical temperatures in simple cubic and 
square lattices, that according to the last formula is 
given by Tc 3 ^ /Tj 2d ^ = 3/2, is approximately fulfilled 
for U > 20 t. 


B. Compressibility and double occupancy 

The electronic conductivity is an important character¬ 
istic of solid materials. However, in systems of ultracold 
atoms in optical lattices, it is easier to measure the com¬ 
pressibility, which is qualitatively related to the conduc¬ 
tivity: an electronic system is conducting if the many- 
body state of the electrons is compressible. This elec¬ 
tronic compressibility is defined as the derivative of the 
electron density with respect to the chemical potential, 
K e = Fig. pphows the comparison of the contour lines 
with K e = 0.017* for particular values of t z . The shift of 
the contour line towards larger interaction strengths for 
larger t z in the unordered phase shows that at certain pa¬ 
rameter values of U and T it is experimentally possible to 
change the compressibility of the system in a wide range 
by just varying the hopping parameter t z (i.e. changing 
the laser intensity in one direction) and keeping U and 
T constant. 

The average on-site double occupancy given as Di = 
is another useful indicator for the degree of lo¬ 
calization of the fermions and an observable quantity of 
great importance in optical lattice experiments 01 El ED . 
In the regime of strong coupling and low temperatures it 


was found to exhibit an increase in the presence of an¬ 
tiferromagnetic correlations that appear with decreasing 
temperature [221122] ■ In order to analyze the relationship 
between the double occupancy and magnetic ordering, we 
plot both D and m as a function of temperature in one 
diagram (Fig. [4]). We clearly see that the onset of stag¬ 
gered magnetization coincides with an enhancement in 
double occupancy at large interaction strengths, in ac¬ 
cordance with the results given in 22 ■ We also find that 
at given system parameters for the two-dimensional case 
(t z /t = 0) there is a good agreement with the results ob¬ 
tained by means of numerical linked cluster expansions 
technique in |23] . In the three-dimensional case, however, 
our observations show stronger dependencies in compar¬ 
ison with more accurate dynamical cluster approxima¬ 
tion results [21] ■ At intermediate coupling, the signal be¬ 
comes less pronounced, and below a certain value of f7, 
we even observe converse behavior: the double occupancy 
smoothly increases just before the system enters the an¬ 
tiferromagnetic phase and then decreases with increasing 
magnetization. Comparing these diagrams for different 
t z we note that the effect of the double occupancy en¬ 
hancement due to antiferromagnetic correlations occurs 
at larger U for larger t z values, e.g. at U ~ lOf for t z = t 
and at U ss 8t for t z = 0. These values correspond to the 
paramagnetic crossover region between the metallic and 
Mott insulating states, which is shifted to lower interac¬ 
tion strengths for smaller t z - values (see also Fig. [3]). 

In general, we observe higher double occupancy at 
lower interaction strength. This effect is directly un¬ 
derstandable from the model Hamiltonian ([!]): double 
occupation of a particular lattice site becomes energeti¬ 
cally unfavorable in case of high interaction “costs”. The 
increase of double occupancy in the antiferromagnetic 
phase at strong coupling is an effect of virtual hopping: 
In the antiferromagnetically ordered system each parti¬ 
cle is surrounded by particles with an opposite spin, al¬ 
lowing it to virtually hop to all Z next neighbors and 
by this to lower its energy, whereas in the paramagnetic 
regime only 50% of the next neighbors have the oppo¬ 
site spin (see also Ref. j22] for details). In the regime 
of weak coupling, the paramagnetic state is metallic, i.e. 
the particles are itinerant. Thus, the double occupancy 
is generally large. However, the temperature affects the 
metallic quality of the state: at higher temperatures, the 
delocalization of particles is disturbed by thermal fluc¬ 
tuations. Therefore, the delocalization and, respectively, 
the double occupancy increases at lower temperature val¬ 
ues. The peculiar cusps in the curve shape emerge at the 
transition to the antiferromagnetic phase: the double oc¬ 
cupancy decreases because the antiferromagnetic order¬ 
ing opens a gap in the charge excitations spectrum, i.e. 
the system becomes an insulator. 

To analyze how the system’s dimensionality affects the 
double occupancy we chose particular pairs of points 
in the magnetic phase diagram, which are indicated in 
Fig- [5] and found a pronounced dependence of D on the 
value of t z . The increase of D is particularly strong when 















5 



Temperature T/t 


FIG. 4. (Color online) Double occupancy (solid lines) and 
magnetization (dashed lines) plotted as a function of temper¬ 
ature for different values of the interaction strength U and 
t z /t = 0, 0.5,1. At large interaction strengths a pronounced 
increase of double occupancy is observed below the Neel tem¬ 
perature. 




the change in t z from zero to t induces a phase transi¬ 
tion, as it is the case in the points 1 (AFM to PM) and 
3 (PM to AFM). In the points 2 and 4-6 the phase does 
not change for any value of t z £ [0, t\, thus the curves 
have nearly the same slopes. The resulting dependence 
of the double occupancy on the hopping parameter t z 
can be used to identify the presence of magnetic cor¬ 
relations or to measure the temperature of the system, 
since D is an experimentally accessible observable, t z is 
a tunable parameter, and the interaction strength U is 
a known quantity of the experimental setup. However, 
it was pointed out in [23] that the increase of double oc¬ 
cupancy in the considered parameter regime occurs not 
only in the Mott-insulating core, but also in areas with 
n < 1. Thus, to identify the increase of D as a signal 
of onsetting antiferromagnetic order in experiments, it 
has to be additionally ensured that the density in the 
measured region of the trap is close to n = 1. 


C. Entropy analysis 

In experiments with ultracold atoms the system is 
characterized in a natural way by the total particle num¬ 
ber N and the total entropy S. These quantities can 
ideally be kept fixed, since the system is isolated from 
the environment and the experiments are assumed to be 
carried out adiabatically. The entropy per site, however, 
is variable and is often used as a thermometer, since it is 
independent of the trap frequency, while the temperature 
changes even when the trap frequency is increased adi¬ 
abatically [25]. To determine the entropy theoretically, 
we use the definition from statistical thermodynamics: 
S = In Q (ks = 1), where is the number of possible mi¬ 
crostates of the system under consideration. In the Hub¬ 
bard model within the single-band approximation 0- a 


FIG. 5. (Color online) Points in the magnetic phase diagram 
with transition lines between the PM and AFM phase for 2d 
(t z /t = 0) and 3d (t z /t = 1) and analysis of the pairs of 
points from adjacent regions. Double occupancy is plotted 
against t z . For each curve, the values for T and U are kept 
constant. (1): U = 2.5t,T = 0.07i, (2): U = 10 t,T = 0.2 1 , 
(3): U = 15t, T = 0.35f, (4): U = 2.5f,T = 0.25f, (5): 
U = 17.5 1, T = 0.4t, (6): U = 101, T = 0.55 1. 


single site can be found in four different states, thus the 
maximal entropy per site is given by s max = ln(4). This 
limit is reached at half filling (the number of particles 
coincides with the number of sites) and vanishing inter¬ 
action U = 0, where all four possible states have the 
same probability in the framework of the microcanonical 
description. In case of strong interactions, however, the 
system becomes a Mott insulator with only two possible 
states (|t) and ||)), such that the upper entropy bound 
is given by s = ln(2) « 0.69. The calculations with the 
combined DMFT+LDA method provide us with expec¬ 
tation values of the local particle density (“filling”) n at 
different values of the chemical potential. Given these 
quantities, we can use a Maxwell relation for the entropy 
per site m 


ds 

<9/r 


dn 

fff 


s{no,T) 




( 8 ) 


Fig. § shows the isentropic curves in the T /[/-plane ob¬ 
tained for different values of the hopping parameter t z . 
Naturally, larger entropy values are found at higher tem¬ 
perature. We observe a characteristic shape of the isen¬ 
tropic curves for small entropy values (s < 0.7), which 
have a negative slope at weak coupling, bend sharply 
when entering the antiferromagnetic region and then 
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FIG. 6. (Color online). Isentropic curves in the T/CI-plane 
obtained by LDA+DMFT at half filling (po = Uj 2) for the 
hopping parameters t z /t = 0, 0.6,1. 


follow the curvature of the magnetic phase transition 
boundary (black curve). On the left side of the phase 
diagram the so-called Pomeranchuk cooling effect can 
be observed: an increase of the interaction strength U 
at constant entropy leads to a decrease in temperature 
beyond non-interacting band structure effects |26| . The 
isentropic curves are shifted towards larger values of T 
and U when t z increases, which makes the region with 
the Pomeranchuk effect larger. However, according to 
an additional analysis performed we can conclude that 
this enlargement is mainly based on the increase of the 
bandwidth. Comparing our results in the limiting cases 
t z /t = 0 and t z /t = 1 to references [23] and [24] (27], re¬ 
spectively, we observe nearly the same quantitative over¬ 
estimates of critical entropies by the DMFT method. 
At the same time, the qualitative behavior persists and 
agrees well in both limits, and therefore we expect that 
the results in the dimensional crossover region remain of 
a high importance. 

To analyze the effects of hopping anisotropy, we plot 
the temperature, scaled by the bandwidth W (in order to 
exclude its broadening effect), as a function of the hop¬ 
ping parameter t z for particular entropy values. Fig. [7] 
shows these plots for four different interaction strengths. 
We observe that, in general, in the 2d-case {t z /t = 0) the 
rescaled temperature at a given entropy value is larger 
than in the 3d-case (t z /t = 1). This effect can be easily 
understood when we consider the entropy as a measure 
for disorder: the transition from 2d into 3d increases the 
number of the nearest neighbors from Z = 4 to Z = 6, 
opening two more tunneling options for a single site (i.e. 
the positive and the negative z-direction). Hence, the 
entropy of the system should increase, but since its value 
is kept fixed, its conjugated thermodynamic variable - 
the temperature - decreases. However, the decrease of 
T/W is not monotone in general: for particular U- and 
s-values, the isentropic curves exhibit a minimum at in¬ 
termediate t^-values. This effect is based on the magnetic 
phase transitions of the system. The slope of the isen¬ 
tropic curves becomes positive when the initially antifer- 
romagnetically ordered system (regions shaded in grey 
in the figure) turns into the paramagnetic state. Again, 
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FIG. 7. (Color online). Temperature T, scaled by the band¬ 
width W = 8t + 4 t z , is plotted as a function of the hopping 
parameter t z /t for fixed values of entropy s at different inter¬ 
action strengths Ujt = 0, 3.5, 6,8.5. 


this effect can be explained in the entropy picture: In 
the antiferromagnetic phase the spins are ordered in a 
checkerboard pattern, so that hopping is possible to each 
of the nearest neighbor sites. In the paramagnetic phase, 
however, on average only 50% of the nearest neighbors 
have an opposite spin. Thus, with the transition into the 
paramagnetic state the entropy should decrease, and by 
the same argument as above the temperature increases. 


D. Signatures of anisotropy-driven transitions in 
the real-space density profiles 

The real space density distribution of the atoms 
is experimentally accessible with in-situ imaging tech¬ 
niques HHini na ns]- Within the local density approx¬ 
imation this is equivalent to measuring the density as 
a function of chemical potential for a given tempera¬ 
ture [25]. We use the data obtained by the LDA+DMFT 
calculations performed for different parameter regimes, 
which are indicated by points in the upper diagram in 
Fig- i The resulting particle density profiles and the 
magnetization profiles in the trap are shown in below. 
Note that at lower temperature the 2d- and the 3d- 
system are in different magnetic phases: In the regime 
of weak coupling U = 3t at T = O.lf the 2d-system is in 
the antiferromagnetic insulating state, whereas the 3d- 
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system is in the region of the paramagnetic metal (al¬ 
though very close to the magnetic transition line, thus a 
finite magnetization is present). Similarly, in the regime 
of strong interaction U = 13 1 at T = 0.37 1 the 2d-case 
corresponds to the paramagnetic insulator, while the 3d- 
case is in the antiferromagnetic insulating phase. In each 
case the particle density equals one in the trap center and 
decreases with increasing distance from the trap center 
to finally vanish at certain values of r. The shape of 
the density distribution differs depending on the param¬ 
eter values U and T. In the regime of strong coupling 
(U/t = 13) both density profiles exhibit a wide plateau 
in the trap center, indicating an insulating state. How¬ 
ever, in case of weak interaction this plateau can be iden¬ 
tified as a signature of the antiferromagnetically ordered 
state: it only appears when significant staggered magne¬ 
tization is present. Thus, in the regime of weak coupling 
the real-space density profile can be used in the experi¬ 
ment to detect an antiferromagnetic phase without mea¬ 
suring local quantities. However, in other regimes it re¬ 
mains essential to access local observables such as magne¬ 
tization and double occupancy in order to determine the 
system’s phase. The sharp kinks, observable in case of 
U/t = 3 ,T/t = 0.1 (2d-curve) and U/t = 13,T/t = 0.37 
(3d-curve), which seem to indicate the appearance of 
the antiferromagnetic phase, are known to be artifacts 
of LDA+DMFT and are smoothed out in the real sys¬ 
tem by proximity effects. 


V. CONCLUSIONS 

We studied the effects of hopping anisotropy on the 
physical properties of fermionic atoms in optical lattices 
with simple cubic geometry. The analysis is based on 
calculations for the hopping parameters t x = t y = t and 
t z /t £ [0.0,1.0] performed with the DMFT approach. We 
found that all characteristic many-body properties and 
quantum phases of the system depend on t z . In particu¬ 
lar, we analyzed the significant changes in the magnetic 
phase diagram that depend on the interaction strength 
in a non-trivial way: the critical temperature for antifer¬ 
romagnetic long-range order decreases with t z at small 
values of U, but increases with t z in the regime of strong 
coupling. This behavior is consistent with the analytical 
expressions from the mean-field approaches for the lim¬ 
iting cases U Zt and U <C Zt, which have different 
dependence on the coordination number Z. 

The analysis of metallic vs. insulating properties 
showed that it is experimentally possible to change the 
compressibility of the system by just varying t z and keep¬ 
ing the other parameters (U and T ) constant. The dou¬ 
ble occupancy increases with increasing t z and shows a 
pronounced dependence on the magnetic phase of the 
system. Thus, in the experiments the double occupancy 
can be used to determine the presence of magnetic cor¬ 
relations in the system or to measure the temperature. 
However, we note that the obtained results for the dou¬ 



FIG. 8. (Color online). A selection of points in the mag¬ 
netic phase diagram for the analysis of the particle density 
distribution in a harmonic trap and its dependence on the 
quantum phase in the system. Below: Particle distribution 
in real space (solid line) and magnetization (dashed line) for 
tg/t = 0 (blue, thin) and t z /t = 1 (red, thick), in the regimes 
of weak (U = 3 1) and strong (U = 13f) coupling at different 
temperatures. On the z-axis is the radial distance from the 
trap center, measured in units of the lattice constant a = 1. 
The trapping potential is V — 0.1 1. 


ble occupancy are based on DMFT, i.e., momentum- 
independent analysis. Thus, the influence of non-local 
correlations close to transitions is not properly accouted 
for and this can lead to some deviations between theo¬ 
retical predictions and experimental results. Therefore, a 
further improvement in this direction would be a theoret¬ 
ical analysis of the double occupancy by more accurate 
approach, e.g. the dynamical cluster approximation m 
or diagrammatic determinant Monte-Carlo method [27], 
in the region of the studied dimensional crossover. 

In combination with the local density approximation 
(DMFT+LDA), we analyzed the entropy and the real- 
space density profiles of the system with a harmonic 
trapping potential. In the analysis of the entropy we 
find that the parameter region with the presence of the 
Pomeranchuk effect becomes larger with increasing t z . 

















































































At constant entropy, the rescaled temperature generally 
decreases with t z , which is simply explained by the in¬ 
crease of the coordination number. However, when the 
initially antiferromagnetically ordered system turns into 
the paramagnetic state, the temperature rises again. Our 
analysis of the particle density distribution in real-space 
showed that at weak coupling the transition into the mag¬ 
netically ordered state corresponds to the appearance of 
a plateau in the density profile. Thus, this signal could 
also be used to detect the antiferromagnetic phase in the 


experiment. 
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Appendix A: Non-interacting density of states for 
anisotropic systems 

The general expression for the density of states (DOS) 
in n dimensions is given by 

D n(e) = J d n k 5{e - e(k)) (Al) 

BZ 

with integration over the Brillouin zone. The energy is a 
function of the wave vector k. For cubic lattice geometry 
its general form is given by 

e(k) = —2 ^ t a cos (k a a a ), 

at 

where t a denotes the hopping parameter and a a the lat¬ 
tice spacing in a-direction, and k a are the corresponding 
components of the wave vector. We consider lattice ge¬ 
ometries with equal lattice spacings in all directions and 
set a a = 1 for simplicity. 

Our aim is to calculate the non-interacting density of 
states for anisotropic hopping, i.e. different values of t a 
for a = x, y, z , in particular, for the case of isotropic hop¬ 
ping in the xy-plane and varying the hopping parameter 
in z-direction. In the following, the closed analytical ex¬ 
pression for the density of states is derived in terms of 
complete elliptic integrals of the first kind. 
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1. DOS in One Dimension 


In one dimension the energy and the density of states 
are given by 

e(k x ) = —2t x cos(k x ) = -e 0 cos(fc x ), 

7 r 

D i( e ) = :^ r j dk x 5 (e - e{k x )). 

— 7T 

The integral is evaluated using the following property of 
the Delta function: 


with g(xi) = 0 . 


(A2) 


Due to the periodicity of the system the integration 
limits in the integral for D\ can be extended to infinity 
and the equation for the DOS reads 




1 


2tt 

1 


— Eo sin(arccos(—e/e 0 ))| 

1 


7 T | — eq sin(arccos(— e/eq))\ 


The sum just gives a factor of 2, since the two roots 
of the function g(k x ) = e + EqCO s(k x ) are given by 
± arccos(— e/eo), in accordance with the symmetry of the 
cosine function. In the last step we use the relation 

sin(arccos(x)) = cos(arcsin(x)) = \/1 — x 2 , (A3) 


and the final result reads: 


Di(e) 


1 

7TV^O - £2 


2. DOS in Two Dimensions 


In two dimensions, we consider a square lattice with 
anisotropic hopping t x ^t y . The energy is given by 

e(k x , k y ) = —2 t x cos (k x ) — 2 t y cos{k y ) 

and the density of states reads 


7T 7T 


( 2 tt) ; 


D 2 (£) — ^27r)^ f dk x J dky 5(e s(^k X: ky^ 

— 7T —TV 

7r 7T 

J dk x J dky 5 (e + 2 t x cos (k x ) + 2 t y cos (k y )) 


—7r —7r 

7r 7r 


2 ^ - J dk x j dky £(£+ Acosfe) + cos(fc y )), 


0 0 


where in the last step we introduce dimensionless param¬ 
eters i = £ j 2t y and A = t x /t y and change the integration 
limits according to the symmetry of the integrand. In the 
next step the substitution a = cos (k a ) is performed for 
a = x,y and the DOS becomes: 


D2{E) 



dy 


8{e + xA + y ) 
y/(l-x*)(l-y*Y 


The evaluation of the (5-function for y according to 
Eq. (|A2| gives 


D 2 (e) = 


1 


7t 2 2L, 


dx 


1 


a /(1 - x 2 ) [1 - (e + x A) 2 ] 


The integration limits of x are changed by this trans¬ 
formation and are now determined by the square root 
in the denominator. Using the relation (A3) we obtain 


y/(l — x 2 )[l — (i + xA) 2 ] = sin(arccos(a;)) sin(arccos(e+ 
xA)). The argument of arccos is restricted to values in 
the interval [- 1 , 1 ], thus the two conditions 


— 1 < x < 1 and — 1 < (e + xA) < 1 


have to be satisfied. This leads to the following limits: 


Xmin = max 


- 1. 


— 1 — £ 

A 


1 — £ 


with —1 — A<£<l + A. Evaluation of the dependence 
on £ leads to four distinct regions: 


£ > 0 

I: 1-A<e<1 + A 

-1 < x < ^ 


II: 0 < £ < 1 - A 

-1 < X < 1 

£ < 0 

III: A - 1 < £ < 0 

-1 < X < 1 


IV: 0 < £ < A - 1 

—< x < 1 





But due to the symmetry of the integrand we find that 
D j (£ > 0;A) = D iv (e < 0; A) and D ti {e > 0; A) = 
D 111 (i < 0; A). This result will be considered after the 
following transformation. 

To evaluate the integral we factorize the polynome 
under the square root in x, bringing it to the form 
p {x) = IL (x — Xj), with Xj being the roots of the poly¬ 
nome, 

1 — £ —1 — £ 

Xi = 1, x 2 = -1, x 3 = —£~, x 4 =———, 

so that the DOS reads: 

If 1 

D 2 (e) = -t—— / dx — , 

21 ' n 2 2t y J y /p(x) 

In the next step, the following substitution is applied 
to the factorized polynome: 


- 6)-5(P- 7 ) sin 2 (<f>) 
{13-5)- (fi-i) sin 2 ( 0 ) 


(A4) 
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with 6 < 7 < j3 < a being the real roots of the polynome. 
This way, the integral is transformed into an elliptic in¬ 
tegral of the first kind: 


dx 


d<j) 


Vm l)(P J \Jl — m sin 2 ((f>) 


with the parameter m given by 

(/3 — 7)(a — 6) 


7 < x < /3. 


~S) 

The integration limits transform according to the formula 


4>(x) = arcsin 


'(x - 7)(/3 - 5) 
(x — 5)(/3 — 7) 


This transformation is performed separately in the two 
regions stated above, i.e. for D 1 = D 11 and D n = D 111 , 
since the grading of the roots Xi is dependent on e. For 
the regions I and IV the grading is given by Xi > X3 > 
X2 > X4 and the DOS transforms into 


Dz(e) = ^k K(m ° 


(1 + A) 2 - e 2 
4A 


and for the regions II and III we obtain x 3 > X\ > x 2 > 
X4 with the DOS given by 


D 2 (e) = 


2A ■ K(m b ) 


tt 2 2 t x y/(\ + A) 2 + e 


m b = 


where K(m) denotes a complete elliptic integral of the 
first kind with the parameter to, given in the Legendre 
form by 


tt/2 

A'(to-) = f —j= 

0 V 1 


dcj) 


— to sin 2 (</>) 


(A5) 


3. DOS in Three Dimensions 

In three dimensions we consider a simple cubic lat¬ 
tice with equal hopping parameters in the aiy-plane, 
t z 7^ t x = t y = t. The kinetic energy is given 
by 

e(k x , k y , k z ) = —2t[cos(k x ) + cos(k v )} - 2 t z cos (k y ) 
and the density of states is 

7T 7T 7T 

— ~(2/7 r)^ f J J 5(s s(kxi ky : k z y) 


— 7T —7T —7T 

7T 7T 7T 


1 


- j dk x J dk y J dk z 
000 


x 5(e + cos (k x ) + cos (k y ) + A cos(fc 2 )), 


where we introduce the dimensionless parameters s = 
e/2 1 and A = t z /t and reduce the integration range using 
the symmetry of the integral. Proceeding in an analogous 
way as in 2d, we first substitute the cosine terms by a = 
cos (k a ) for a = x, y, z to obtain 
111 

1 f , f . f . S(e + x + y + zA) 


D 3 (e) = 


dx dy dz 


n 3 2tj J yj{l-x 2 )(l-y 2 )(l- z 2 ) 

and evaluate the 5-function for y according to Eq. ( |A2| ), 
taking x and z as constants, which leads to 

Ds(e) = 

f dx [ dz 1 

^ 3 2 1 J J ^/(l - x 2 )(l - z 2 )[l - (e + x + zA) 2 ] 

Expressing the denominator in terms of the rela¬ 


tion (A3) and considering the domain of arccos(a;) as 


shown in the previous section, we find the following inte¬ 
gration limits: 

1 / lj Xmin / X ^ X max , 

x min = max[—1, -1 - (e + zA)], 

Xmax = min[l, 1 - (s + zA)}. 

The dependence on (e + zA) leads to two distinct regions 
that need to be considered separately in the following 
discussion: 


4A 

(e + zA) > 0 

I: — 1 < x < 1 — (i + zA) 

(1 + A) 2 - e 2 ’ 

(s + zA) < 0 

II: 

— ! — (£ + zA) < x < 1 


Next, the polynome under the square root is factorized 
in x. The roots of the polynome function are 

X\ = 1, x 2 = — 1, X 3 = — 1 — e — zA, £4 = 1 — e — zA, 

and the DOS becomes 

1 /',/■, 1 


D 3 (s) = 


TT 3 2t 


dx / dz 


'(i-*2)n(*-*o 

i =1 


Now we perform a substitution given in Eq. (A4) to 

transform the DOS into 

n , X 1 f dz 2 

D-iie) = —— / — . — . x 

n 2t J z<2 ) V( a ~ t)(/3 - 6) 

f d(f> 


\Jl — to sin 2 (</>) 


Here, the two regions shown above need to be considered 
separately. In the region I the grading of the roots reads 
x\ > X4 > x 2 > X3 and in region II we find X4 > aq > 
x 3 > x 2 , but both lead to the same transformation, thus 
the DOS is given by 
1 


D 3 {e) = 


dz 


7T 3 2 1 J y/l _ 

-1 


-K(m), to = 1 — 


£ + zA 


with the compete elliptic integral AT(to) defined in 
Eq. (1A5). 
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Appendix B: The AGM method 

Complete elliptic integrals of the first kind can be cal¬ 
culated exactly by the Gaussian AGM method. 

It was found by Gauss that the sequences for geo¬ 
metric and arithmetic mean have the same limit - the 
arithmetic-geometric mean. Given two starting values 
clq and b 0 , the sequences read: 

ai = -(ao + bo ), ..., cin+i = -(d/v ^n) 

and 

bi = \Jaob 0 , ..., 6 jv+i = \J• ajsibN■ 

In the limit N — > oo both sequences converge to the same 
value: 

lim aw = lim bjy = agm(a,b). 

N—y oo N—>oo 


For the integration limits we solve the quadratic equa¬ 
tion for t 


ti /2 = x ± \J x 2 + ab , 


which gives two equivalent results, since the integrand 
function is symmetric. Choosing t\ we obtain 


x —> oo : t\ —► oo, x -A —oo : ti —> 0, 
thus the final result reads: 


/ dt 

/,.o . 2 w, o_ s _ h 2\ = J ( aAr ’ bN ) 


For elliptic integrals of the first kind, Gauss derived 
the simple relation 


This yields the general relation 


OO 

I(a,b) = J 


dx 


sj (x 2 + a 2 )(x 2 + b 2 ) 2 • agm(a, b )' 


I(dN+l, bN+l) = I ^2 ( aN + ^n) , \f(lNbN 


To show that, we consider the integral 


I( a N+i,b N+1 ) — - [ —j= 

2 -L /( 


dx 


(x 2 + a 2 N+1 ){x 2 + b 2 N+1 ) 


and perform the substitution 


1 , 

x = - \t — 


ClNbN 


dx = - I 1 + 


t ' 2 J 


dt. 


The integral then reads 


From this, it directly follows that /(aAq&jv) is indepen¬ 
dent of N: 


I(a 0 ,b 0 ) = I(a N ,b N ) = lim I(a n ,b n ) = I(c,c) 

n—> oo 


with c = agm(a, b). Note that J(c, c) is directly solvable: 


I{o,n+i, bu+i) 

-u— 


(1 + a^bN/t 2 ) dt 


(/ - +a 


at+i 


+b% + 1 


OO 

/(c,c) = J 


dt 


dt 


yj(t 2 +c 2 )(t 2 + c 2 ) J t 2 +c 2 


7 r 

2c' 


— /" cZt_ 1 + a^bx/t 2 _ 

^ ^[®2 (f 2 + a^)(f 2 + 6^)] [i< 2 (1 + aArkjv/f 2 )] 

r _dt_ 

J yj ( t 2 + a 2 N ){t 2 + b 2 N ) 


The relation between the Gauss form for the complete 
elliptic integral of the first kind and the Legendre form 
K(m) can be found by the substitution tan cf> = x/b. 



